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An optimal shape design problem for plates 


Klaus Deckelnick* Michael Hinze"^ Tobias Jordan^" 


Abstract 

We consider an optimal shape design problem for the plate equation, where the variable 
thickness of the plate is the design function. This problem can be formulated as a control 
in the coefficient PDE-constrained optimal control problem with additional control and 
state constraints. The state constraints are treated with a Moreau-Yosida regularization 
of a dual problem. Variational discretization is employed for discrete approximation of 
the optimal control problem. For discretization of the state in the mixed formulation we 
compare the standard continuous piecewise linear ansatz with a piecewise constant one 
based on the lowest-order Raviart-Thomas mixed finite element. We derive bounds for the 
discretization and regularization errors and also address the coupling of the regularization 
parameter and finite element grid size. The numerical solution of the optimal control 
problem is realized with a semismooth Newton algorithm. Numerical examples show the 
performance of the method. 

Key words, elliptic optimal control problem, optimal shape design, pointwise state con¬ 
straints, Moreau-Yosida regularization, error estimates. 
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1 Introduction 

This work is devoted to the numerical analysis and solution of an optimal control problem for 
a plate with variable thickness. The state equation 

A{u^Ay) — f in 12 

can be used to model the relation between the (small) deflection y and the thickness u of 
a (thin) plate under the force of a transverse load /. The domain 12 C represents the 
unloaded plate’s midplane and we assume its boundary to be simply supported, i.e., 

y — Ay = 0 on 512. 

Invoking a pointwise lower bound on the state y and pointwise almost everywhere box con¬ 
straints on the control and minimizing the volume of the plate given by the cost functional 

/ u{x) dx 

Jn 

*Institut fiir Analysis und Numerik, Otto-von-Guericke-Universitffi Magdeburg, Universitatsplatz 2, 39106 
Magdeburg, Germany (klaus.deckelnick@ovgu.de). 

^Schwerpunkt Optimierung und Approximation, University Hamburg, Bundesstrake 55, 20146 Hamburg, 
Germany (michael.hinze@uni-hamburg.de; tobias.jordan@uni-hamburg.de). 


1 



lead to a control in coefficients problem, which can also be viewed as an optimal shape design 
problem. In [Sprekels and Tiba 1999| this optimization problem is analyzed using a transfor¬ 
mation to a dual problem. 

Building upon this duality, it is our aim to solve the control problem with a finite element 
approximation that is suitable with regard to the necessary optimality conditions. To this end 
we compare variational discretization of the control problem (cf. [Hinze 2005| ) based on either 
the lowest-order Raviart-Thomas mixed finite element or piecewise linear continuous finite 
elements for the discretization of the Poisson equation. The pointwise state constraints, which 
are responsible for the low regularity of the Lagrange multiplier, are treated with the help of 
Moreau-Yosida regularization (cf. [Hintermiiller and Knnisch 2006j ). The numerical solution 
to the control problem is computed via a path-following algorithm that simultaneously refines 
the mesh and follows the homotopy generated by the regularization parameter. The resulting 
subproblems are solved by a semismooth Newton method. 

To the best of the authors’ knowledge this is the first contribution to numerical analysis of a 
“control in the coefficients” problem for biharmonic equations including state constraints. The 
mathematical techniques applied in the numerical analysis of the regularized control problem 
are related to the relaxation of state constraints as proposed in IHintermiiller and Hinze 2009] 
and to |Deckelnick, Gunther and Hinze 2009| , where the Raviart-Thomas mixed finite element 
was employed in the context of gradient constraints. 

The present work is organized as follows. In Section the optimal control problem and 
its dual problem are introduced. The regularization of the dual problem is investigated in 
Section Section deals with the discretization of the regularized problems and with the 
related error bounds. Finally, in Section the original control problem is solved with a 
Newton-type path-following method. Numerical examples are presented which validate our 
analytical findings. 

2 The optimization problems 

Let C d G {2,3}, be a bounded domain with smooth boundary 90. The Dirichlet 
problem for the Poisson equation 

-Ay = g in O, 
y — 0 on dVt 

admits for every g G a unique solution y \—T{g) H satisfying 

\\y\\H^{Q) — ^ II^IIl2(o) • (2-2) 

In order to define the control problems considered in this paper we introduce the admissible 
sets for controls and states according to 

f^ad •= {u ^ \ m < u < M a.e. in } , 

^ad := { / e L^{n) I M-^ < I < m-^ a.e. in LJ } , 

Y^d — { 2/ ^ C{n) \ y>-r in } , 

where r > 0 and 0 < m < M are positive real constants. For a given / G L?{^) we consider 
the following optimal control problems (cf. [Sprekels and Tiba 1999| , problems Pi and Di): 


( 2 . 1 ) 















subject to 


A{u^Ay) = f 

in D, 

(2.3) 

y = 0 

on 5D, 

(2.4) 

1 

> 

II 

0 

on 5D, 

(2.5) 


y ^ 

denoted as the primal problem (P), representing the physical control problem motivated in 
the introduction, and secondly, with the datum z E V induced by / 


—A 2 ; — f in fi, 

^ = 0 on 

the dual problem (D) 

min Jil):— [ dx 

ieLo-(Q) 

subject to 

—Ay — zl in fi, 

y — d on 

^ ^ -^ad? 

y ^ ^d? 


(D) 


( 2 . 6 ) 

(2.7) 


which is analytically and numerically advantageous, in that it is convex and contains two 
coupled second order equations instead of the fourth order equation (2.3). It will therefore 


serve as a basis for our analysis in the remaining sections. For every u E ?7ad fhe system (2.3)- 

V with v?Ay E V. Due to 2 : G L^(D) we 
G P to the dual system (|2.6|)-(|2.7|). 


(2.5) has a unique weak solution y — y{u) E 
have zl E L^(D), and there is a strong solution y = y{l) 
We impose the following Slater condition: 


^Us E C/ad,^5 >0: ys = y{us) > -r + £s, (2.8) 

and recall from [Sprekels and Tiba 1999| that for each pair {y^u) admissible for (P), the 
pair (y,/ = u~^) is admissible for (D) with the same cost and vice versa. Moreover, there 
exists a unique solution ?iopt of problem (P), and /opt = '^op^t unique solution of (D). We 

denote the associated state by ^opt- 

Next we derive optimality conditions characterizing /opt- For this purpose we intro¬ 
duce ^(D), the space of regular Borel measures, which equipped with norm 

IMU(n) = sup /_/d(-) 

feC(d),\f\<iJn 

is the dual space of (7(D). Arguments similar to those used in |Casas 19861 Theorem 2] yield 
the 
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Theorem 2.1. Let the assumption (2.8) hold. A eontrol I G with assoeiated state y — 

y{l) is optimal for the dual problem (D) if and only if there exist q G L‘^{LL) and v G 
sueh that 



— / qAwdx= / wdiy(x) VwgV, 

Jn Jn 


/-4/3^ (A: — /) dx > 0 V A; G Lad, 

(2.9) 

[ {w -y) dz/(x) <0 Viz; G Yad, 

Jn 

(2.10) 

^ ^ Lad, y ^ Y^d 



are satisfied. 

The variational inequality ( |2.9[ ) can be written as a projection formula 




-3/4 


a.e. in 12, 


( 2 . 11 ) 


wher e P\a^ b] is the orthogonal projection onto the real interval [a, 6]. For later use, we note 
that (2.10) is equivalent to 


z/ < 0, 


{y + r) diy{x) — 0. 


( 2 . 12 ) 


It follows from (2.12) that the support of the measure z/ is concentrated in the state-active 
set { X G 12 I = —T }. In particular z/|^^ = 0, see [Casas 1986] . Furthermore, from 

Theorems 4 and 5 in [Casas 1986] we deduce q G 1Fq^’^( 12) for all 1 < s < djid — 1) 
and q G \ {^ = 


3 Moreau-Yosida regularized problem 


To relax the state constraints, we introduce the Moreau-Yosida regularization (D^) of prob¬ 
lem (D) for a parameter 7 > 0. It reads 

minJ^(/):= [ l-'^/^dx + - [ ((y + r)-)^dx (D^') 

/sLad Jn 2 ^ 

subject to 


y = T{zl). 


Here we set (•)“ := min{0, (•)}. It admits a unique solution P of (D^) with associated state 
denoted by y^ — T{zP). Furthermore, there exists a unique q^ G which together 

with y^ and P satisfies 


/ q^Awdx^ / 7(2 /^ + t) wdx 

JQ Jq 


j (r)-^/^) {k-v) dx>0 


Vw G V 


V fc G LaAt 


(3.1) 

(3.2) 


V G L^A- 
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The term := 7 ( 2 /^ + t) can be regarded as a regularized version of the Lagrange multi¬ 
plier jy G in Theorem 


2.1 


There holds = T{iy^) G V and 


{y^ -w, >0 V w; G Fad- 

This inequality for w G Tad can be argued as follows: 

{y^ -w, u'^)c(n),M(n) = {y^ -w, ^{y^ + T-r-w, 


= (V + ^ , 7 (V + r) + (-(«^ + r ), 

= 7||(V + d“||L2(j^) + (-(«^ + d> >0. 


Now we want to show convergence of the parameterized subproblems (D^) towards the 
unregularized problem (D). We begin with uniform boundedness of primal and dual variables 
with respect to 7 . While the former is obtained immediately through the control constraints 
and ( 2 . 2 ), the latter can be shown as follows. 

Lemma 3.1. Let 7 > 0 and P G Lad be the solution to the problem (D^) with assoeiated 
state and multipliers aeeording to the optimality eonditions. Then there exists a 

eonstant C > 0 independent of 7 sueh that 

ll^^^llLi(n)> < C. 

Proof. To uniformly bound in L^(n) we test ( |3.2[ ) with the Slater element Ig. With the 
help of the adjoint equation (3.1) we get 

c > , V - > ir, z{P - - ys)mn) > 

with a constant C independent of 7 . The desired estimate now follows from 

, V - ys)L^(a) = , y'^ + T-T - ys)L2^a) 

= 7 II (V + r)- ||5^2(j^) + , ys + T)i2(o) 

> 0 + { — p'^ , £s)i2(f2) 


= e,s \ P 




With this bound we want to prove the one on the dual state Let w solve 

—Aw = q^ in fl, 
w = 0 on 


Using (3.1), the embedding of C{Lt) into and the continuous dependence of 


w on q^^ 


we have that 


9^llL2(n) = [ q'^i-^w)dx= j P^wdx <C\\p^\\^i^^^ 

JO JO 

^ C lk"’'|lLl(Q) < C ||(Z^||j;^ 2 (Q) , 




hence ||g^||i 2 (Q) < C. 


□ 
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Next we need to estimate the violation of the state constraint measured in the maximum 
norm with the help of techniques developed in [Hintermuller, Schiela and Wollner 2014| . 

Lemma 3.2. Let P be the solution of problem (D^); the eorresponding state. Then for d G 
{2, 3} we have for every e > 0 a eonstant > 0; independent of sueh that 

Proof We show that Corollary 2.6 of [Hintermiiller, Schiela and Wollner 2014| is applicable. 
To begin with, we note that P = is uniformly bounded in for every s G 

[l^d/{d — 1)), since is uniformly bounded in This implies that zP G is 

uniformly bounded in 7 and that y^ G is uniformly bounded w.r.t. 7 , which by Sobolev 

imbedding theorems holds also in ioi ^ — A — d — e and all 6 : > 0. Thus Corollary 2.6 

in [Hintermuller, Schiela and Wollner 2014| is applicable and delivers our desired bound. □ 

We are now in position to estimate the regularization error. 

Theorem 3.3. Let I and P be the solutions to (D) and resp., with eorresponding states y 

and y^. Then for every 6: > 0 there exists a eonstant > 0^ independent of for whieh it 
holds that 


^^IIl2(0) + 11?^ + II 2 / 2(^ 4 )+^. 


Proof Using I as test function in (|3.2|) and P as test function in (2.9) we obtain 


C\\l-P\\l2^n) 

< {q^ -q, zl -zP)^ 2 ^Q) = - {q^ -q, A{y - 
= (y - V , 7(V + r )-)- (y - V , ’^)c{n),^{n) 

= {y + r, 7(V + r )-)- (V + r , 7(V + r)“) ^2^^) + iv'" - V ^ ^)c{n),J^{n) 
=: {I) + {II) + {III). 


Since y € Y^ad we have (I) < 0. Moreover {II) = —7 ||(y''' + r) ||^ 2 (q) < 0. The third addend 
is treated with the complementarity condition (| 2 . 12 |) for the multiplier v: 


{III) = {y'r + T, + {-y-r, z^)c(Q),^(Q) 

= (y'’' + r , Oc(n),^(Q) =0 

< ((y'^ + r) , Pc{h),.4^{h) • 

With the help of Lemma |3.2| we arrive at 

C\\l-P\\l2^a) < + ^>C(Q),^(Q) ^ l|(y^ + ^)”llL-(Q) ll^lU(n) 

< C’, ( 7 ^/ 4 -!+®) , 

with 7 -independent constants C, (7^. The continuous dependence of the states on the controls 

and the continuous embedding ^ C{f}) allow to extend this estimate to \\y — 

and lly-yT'll^oo(a). □ 
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4 Finite element discretization 


4.1 


Mixed piecewise constant versus piecewise linear approximation 


In this section we turn to the variational discretization of the regularized control problems, 
taking into account the structure imposed by the optimality systems, especially the projection 
formula (2.11) and its discrete counterparts (4.7) and (4.15). The function (P[^4 
applied to the product of two state variables is evaluated with little effort if those variables 
are approximated piecewise constant and yields an implicit piecewise constant discretization 
of the optimal control. Further, the finite element system of the semismooth Newton method 
in Section]^ in particular the parts (A.2) and (A.3) involving the projection formula and its 
generalized derivative, in this situation is easily assembled exactly. 

Approximating the states with piecewise linear, continuous finite elements delivers a more 
involved variational discretization of the controls, since the projection formula then no longer 
implies a piecewise polynomial discretization of the control variable, but rather the negative 
power of the pointwise projection of a piecewise quadratic function. Moreover, the approximate 


computation of the terms (A.4) and (A.5) introduces an additional error. On the other hand, a 


piecewise linear ansatz delivers the higher approximation order two for the states, as opposed 
to an order of at most one for a piecewise constant ansatz. This is supported by the convergence 
rates w.r.t. the grid size h in the error plots in Section]^ and also allows for a better resolution 
of the control active sets. 

In the remainder of this section we give estimates for the overall error in both discrete 
approaches. 


4.2 Variational discretization of (D^) with mixed finite elements 


Following the above remarks, we use a mixed finite element method based on the lowest-order 
Raviart-Thomas element. To begin with we recall the mixed formulation of the Dirichlet- 
problem for the Poisson equation, i.e., for g G y = T{g) and v = Vy there holds 


/ V • w dx + / y div w dx = 0 

Jq Jq 

/ (j) div vdx+ / g (l)dx = 0 

Jn Jn 


Vw G i7(div, fi), 
ycpeL^iCi), 


(4.1) 


where H{div,Q) := { w G | div w G }. For a given right-hand side g G 

we represent the solution of this mixed problem by G{g) := (^, v). In particular, with := V 2 ; 
this means G{f) — ( 2 ;,v^). 

Let a triangulation 7^ of ^ bo given, where h := max^eT/, diam(T) and ^ be the union of 
the elements of with boundary elements allowed to have one curved face. We additionally 
assume that the triangulation is quasi-uniform, i.e., there exists a constant p > 0, independent 
of /i, such that each T G is contained in a ball of radius p~^h and contains a ball of radius ph. 
To define the discrete version of S let us introduce the spaces 


Vfe 

RTo{T) 

Yh 


= RTo{n,7h) := { G H{dw,n) \ ^^h\T e RToiT) 
3 a G 3 /3 G K.: w(a;) = a + jdx 




= I w: 1 

= { G L\n) I VT G Tfe 3/3r G E : <f>h\T = Pt } 


yTe7h}, 

VxgE'^}, 
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For a given g G L^(ri) we set Gh{g) ■= {Vh^^h) ^YhX V/i to be the solution of 

/ Vfe-Wfedx+ / div Wfe dx = 0 Vw/iGVfe, 

./a ./n 

/ 4>h div V/, da; + / gr 0/j da; = 0 
The resulting error satisfies (see [Brezzi and Fortin 1991] ) 

l|y “ yh\\L 2 (a) + l|v - v/i||^ 2 (fj)ci < Ch (||y||Hi(a) + l|v|iiji(a)ci^ 

^ ||2/||/f2(^) < Ch ||^||^2(^) , 

as well as, if ^ G the pointwise estimate (see [Gastaldi and Nochetto 1989[ Cor. 5.5]) 

\\y-yh\\L^{n) + l|v-v/,||^oo(^)d < Ch\\ogh\ ll^ll^oo(^). (4.3) 

The load / induces a discrete datum Zh G via {zh^^z.h) — Gh{f)- 




(4.2) 


Remark 4.1. For our error analysis we require that is bounded uniformly in h in 

both of the considered discretization approaches. 

This is satisfied, e.g., if / G L^(i4), since then ||^ — 


0 by (4.3) and (4.11), resp. 


Note that this regularity restriction is not essentially necessary, cf. [Deckelnick and Hinze 20T4| 
Lemma 3.4], which holds analogously for the scalar states in both discrete approaches and 
allows for / G L^(fl), p > 2, with smaller powers of h in (4.3) and (4.11). 

Let Tad,/i ^ Yh \ (ph\T ^ yT G T/i }. The variational discretization (D^) of 

the regularized control problems (D^) reads 




(»;) 


subject to 


(yh^^h) = Gh{zhl), 

I ^ ^ad- 


We note that (ID^) is still an infinite-dimensional optimization problem similar to (D^) since 
the control I is not discretized. It admits a unique solution /^, which is characterized by the 
optimality system 


/ 

Jq 


(Vh^^h) = Gh{zh ll), 

iyl^lh) = Gh ( 7 (yZ + ^)”)> 

ll e i^ad- 


G -^ad) 


(4.4) 

(4.5) 

(4.6) 


Condition (4.6) is equivalent to the projection formula 

ll = (-P[m4,M4] {^Olzh)) ^ • 

We denote vl := 7 (y(( + t)“ and similarly to the proof of Lemma 
uniformly in h and 7 : 


3.1 


(4.7) 

one obtains boundedness 


















Lemma 4.2. Let 7 > 0 and G Lad be the solution to the problem (ID^) with state ( 2 /^, v^) = 
Gh{zhl]J- Then there exists an ho > 0 and a eonstant (7 > 0 independent of j and of h sueh 
that 


h\\L^{n) 


< C VO < /i < /iQ, V 7 > 0. 


Proof. Let be the adjoint state and Is the Slater element with corresponding discrete 

state {ys^h: ^s,h) — G^izhls), which is a discrete Slater state for all 0 < /i < /iq with some ho > 0 
small enough. In fact, with yg and from (2.8) we obtain from || 2 /s — ^ 0 

that ys^h > + ^s /2 for 0 < /i < ho. 

We test (4.6) with Ig and with the help of the adjoint equation (4.5) and the definition 
of Gh we get 

^ - f 3 ~ 2 ~ ~ ^ ~ ’ 

\ / L^ (^) 

with a constant C independent of 7 and h. We then have 


(7 ’ Vh - ys,h) i2(o) = (7 > yZ + ^ ^ ^2(0) 

= 7 II (7 + r)” ||i 2 (o) + (-7 > 2/3,ft + r) 

>o+(- 7,7) vo</i</io 

V ^ 2/l2(q) 

= fll7Li(a) V0</i</io. 

□ 


Aiming for an estimate of the overall error induced by regularization and discretization we 
apply the approach from [Hintermuller and Hinze 2009] to our problem setting and derive a 
similar asymptotic / 1 , 7 -dependent bound on (/ — Ijf) in which further allows to couple 

the regularization parameter efficiently to the grid size parameter. To this end we need to 
estimate the discretization error for the regularized problems. 

Theorem 4.3. Let P and P^ he the solutions of (D^) and (D^); resp., with eorresponding 
states y^ and ( 2 /^, V/J)- Then there is an /iq > 0 and a 7 - and h-independent eonstant C, sueh 
that for all h G (0, ho) and a// 7 > 0 

II'’ - '2lli-(n) + llw’ - sZlli-M S ^ 

Proof. We define the auxiliary variable (g^,v^) = Gh{i^^) and test the problems’ variational 
inequalities with the respective solutions to obtain 

= {q^z - q^Zh , ll - P) ^2(0) + {q^Zh -q'ilzhJl- P) ^2(o) + (^ft^ft -qlzu^- P) ^2(n) 

=: (I) + {II) + {III). 

In view of Lemma |3.1| we have for sufficiently small 0 < h < ho that 

(I) < Il9’lli=(n, 11^ - 2 /.lli 2 (ni ll'Z - < Ch K - < Ch. 
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For the second addend we find 


m < P/.llL-(n) Il9^ - 9jllL.(n) ll'k - ' 1 lL»(n) £ C Ik’ - «||y ^ 

To estimate the finite element error we set p = sgn(g^ — q^) G L‘^{Q) (cf. the 

proof of Lemma 8.3.11 on p. 228 in [Brenner and Scott 2008| ) and consider 

- Qh ^ p)l^{q) = (q^ ^ p)lHQ) - (qh ^ p)l^q) • (4-8) 

Defining (?/^, v^) = G{p) we have (cf. [Casas 1985| . proof of Theorem 3 with piecewise linear 
finite elements) 


{q^ 5 p)l‘^(q) — {q^ 5 ^y^)L‘^{Q) — ( ^q^ ? y^)L‘^{Q) — (y^ ? 

and furthermore, setting ( 2 /^,v^) = Gh{p) and using the definitions of and 

(Qh , P)LHa) = -{€, = (v^ • , l)^ 2 (o) = “ ( 2 /^ divv);)^ 2 (o) 

Inserting these into ( |4.8| ), we conclude with ||p||/^oo(^) < 1 that 

Ik^ - Qhh^ici) = {(P - Ph > p)l^q) = iP - Ph > sgn(Q^ - gD)L 2 (Q) = [y^ - yl , 

^ \\y^ - yh\\Lo-(Q) ^ Ch\logh\ blkoo(n) < Ch\logh\. 

Finally, we set v/^) = Gh{zh G) and rewrite the third addend {HI) using the definitions 
of q^ and q^ together with the monotonicity of the min-function 

(I/I) = (7(y^ + r)“ - -/{yl + t)“ , yl - yh)^^/^) 

= (7(V + r)- - -/{yl + t)- , yl - y^)+ {-/{y'^ + r)" - -/{yl + r)” , - %) 

-V-^ 

<0 

< (7(V + P~ - livl + r)~ , V - Vh)■ 

Now let (y,v) = G{zhP) and similar to the proof of [Hintermiiller and Hinze 2009[ Theo¬ 
rem 3.5] one has for 0 < /i < /iq 


(7(V + ^) - l{yl + r) , V - yh)L 2 /a) < max 11 
< C’liv -y|lL-(n) + <^l|y-y/illL-(n) < C \\y^ - y\\H2/a) + Ch\logh\ \\zh\\L^/a) k^llL-(a) 




/i. llLi(r 2 ) 


} \\y 


yh\\L°°(n) 


<c\ 


IL°° 


(O) Ik - + Ch |log h\ <Ch + Ch |log h\ < Ch |log h\ 


Altogether, we obtain the proposed err or bo und for the controls. The estimate for the states 
follows similarly to the proof of Theoremjs^and with (4.2) and (4.3). With {p,P) = G{zhll) 
there holds for sufficiently small h > 0 that 


\y'^ yhWi^/n) 


^ IIV y\\L°°{fi) + 


y-y^ 

+ 

s" - vl 


L°°(a) 



y-y^ 


L°°(a) 

+ Ch |logh| |k/,||ioo(n) ll^ftlLoo(n) 


^C'llV y|l_ff 2 (f 2 ) + C y y H^(Q) “ V--, V*. 

< c ||n||^oo(Q) \\z - zh\\L 2 /u) + c lk/,||ioo(Q) ||r - il\\^ 2 /a) + ch \iogh\ 

< Ch ||/||^ 2 (J,) + Cp/^ \logh\C^ + Ch |logh| 

< ChC^ lloghG^ 


□ 
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Combination of the estimates for the two error components immediately gives a bound for 
the overall error. 


Theorem 4.4. Let I and denote the solutions of (D) and (ID^); resp., with eorresponding 
states ( 2 /,v) and Then there exist an ho > 0 and for every e > 0 a and h- 

independent eonstant Cs, sueh that for all 0 < h < ho and 7 > 0 




L2(0) 


+ 


\y-yh\ 


L^{n) 


<a 


^7 2(^ 4)+^ + /2,2 |log/l|2^ . ( 4 . 9 ) 


Proof We split the overall control error into the sum of ||/ — and ||n — 

apply Theorems |3.3| and |4.3| to estimate the regularization and discretization errors. For the 

states we proceed analogously. □ 


Estimate (4.9) now suggests a coupling of 7 and the grid size h of the form 7 = 0{h 
Equilibrating the errors depending on the dimension d, we obtain 



< Ce/12 ^ 


for = 2 , if d = 2 , 
for = 4, if d = 3. 


4.3 Variational discretization of (D^) with continuous piecewise linear finite 
elements 

Next we consider piecewise linear and continuous finite element approximations of the state 
in problem (D^). For g G L?{Pt) given we denote by — Gh{g) E the solution to 

(Vy/,, = {g, ^Wh e U, 

where 


Yh := {Whe C{Cl) I w^T e P\T) VT G Tfe } n Hpn). 


Then, with y — T{g) we have the well known error estimate 


lid d/l|lL 2 ( 0 ) + ^ ll^(d yh)\\L‘^{Q)d < Ch^ Ild|lif 2 (^) < Ch? ||d|lL2(0) • (4-10) 


Provided that g G one can use the L^(ri)-estimate from [Schatz 1998[ Theorem 2.2 

and the subsequent Remark] together with careful control of the constant in the a-priori 


estimate for ( 2 . 1 ) to prove the following bound on the error in the maximum norm 


\\y - yhhoo^a) < Ch^ \loghf ||y||ico(n) 


(4.11) 


We refer to |Deckelnick and Hinze 20081 Lemma 1] for the complete argument. 

Enforcing the state constraints in the (inner) nodes of tho grid, i.e., using the 

space Yad,Ai •= { 0Ai E Y/i I > — T, 1 < < V } of admissible states, the variational 

discretization of the regularized optimization problem (D^) is now straightforward. In prob¬ 
lem (ID^) we only have to replace the discrete solution operator by its counterpart of the 
present solution. With this and the implicit datum Zh = Gh{f) we obtain 


r ■/* +1 (ta+ 


dx 


(Kli) 
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subject to 


Uh — 0 5 

I ^ ^ad- 

Problem (D^ admits a unique solution which together with the discrete adjoint state 
satisfies 


/ 

Jn 


yl = Gh{zhll), 

= Gh {livl + r)“) , 

Ih ^ -^ad- 




(4.12) 

(4.13) 

(4.14) 


The variational inequality (4.14) can again be rewritten as a projection formula 


lh= iPlm\M^]{^qh^h)) 


(4.15) 


We note that {qhZh)\T is a quadratic function, whose projection in general can not be rep¬ 
resented by a polynomial over T. As in Lemma 


4.2 


we infer pi := j{yl + r) < C 

independently of 7 and h. Moreover, Theorem |4. 3 [ holds accordingly. 

Theorem 4.5. Let P and P^ he the solutions of and resp., with eorresponding 

states and y^. Then there exists an /iq > 0 and a 7 - and h-independent positive eonstant C, 
sueh that for all h G (0, /iq) and a// 7 > 0 we have 

11^^ - P\L\u) + l|y^ - ylWHHii) + l|y^ - ylL^ia) ^ • 


Proof. With the adjoint states q^ and qj^ from (3.1) and (4.13), resp., and q^ = we 


obtain as in the proof of Theorem 4.3 
|2 


= {q^z - q^Zh , ll - P) + {q^Zh -q^ZhJl- P) + {q^Zh -qIzhJI- P) ^2(n) 

=: (I) + {II) + {III), 


where, using Lemma 3.1, (4.10), (4.11) and y — T{zhP), yh = Th{zhP) (cf. proof of Theo¬ 
rem |T^, we now can estimate the three addends as follows 


(I) < \\r\Pin) Ik - ZHlPia) T - < Ch\ 

{II) < lk/.lkoo(o) Ik" - qlhnn) P - ^ ’ 

{III) < C Ik - ZhP^^Q) + C' Ik - yfellL-(a) < C {h^ ||/|Il 2 (j^) + |log/i|^ \\zhP\\L 2 (n)) 
< Ch^\loghf, 


with (///) and the estimate of deduced similarly as in the proof of Theorem 

Combining the above estimates completes the proof for the control error. The bounds on' 


state errors can be obtained using similar arguments as for Theorem 4.3 


4.3 

Lhe 

□ 
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The estimate for the overall error follows from combining Theorems 3.3 and |4.5[ 

Theorem 4.6. Let I and be the solutions of (D) and resp., with eorresponding 

states y and y^. Then there exist an ho > 0 and for every e > 0 a 'y- and h-independent 
eonstant C^, sueh that for all h G (0, ho) and a// 7 > 0 we have 

11^ - ^hWmQ) + l|y - yZLi(n) + \\y - ylL^in) < o (7“5(^-2)+" + h |iog/i|). 


Coupling 7 and h again in the form ^ — 0{h to halanee the error eontributions, we obtain 


11 ^ “ ^h\\L^{a) + 



yfiUL^in) 


< Ceh^~^ 


for = 4, 
for = 8, 


ifd = 2, 
if d = 3. 


It turns out that variational discretization with piecewise linear, continuous finite ele¬ 
ments yields a better approximation order for our optimal control problem than variational 
discretization with lowest-order Raviart-Thomas mixed finite elements. In order to achieve 
this, however, a more progressive coupling of regularization parameter and grid size is nec¬ 
essary in the case of piecewise linear, continuous finite elements. While this would mean a 
better approximation of the unregularized solution, it can be a drawback if at the same time 
the problems become harder to solve, cf. the numerical examples in the following section. 


5 Numerical examples 


To approximate the solution of problem (D) we apply a path-following algorithm in the param¬ 
eter 7 , which is sent to oc (cf., e.g., [Hintermiiller and Knnisch 2006j ). The subproblems (ID^) 
and (ID^x)’ solved with a semismooth Newton method, where the parameter 7 

and |4.6[ respectively. The semismooth Newton 
method is described in the Appendix. We conclude this work with two example problems to 
supplement our numerical analysis. 

Example 1: In order to construct an example problem with known solution (cf. Section 2.9 
in [Troltzsch 2010] ). we add the term 

{a/2)\\T{zl)-yn\\l2^^) 

to the cost functional, where a > 0 and yQ G and change the state equation to 


is coupled to h according to Theorems 4.4 


-Ay ^ zl + en, 

with a suitable function oq G see below. To construct the exact solution we set Lt = 

(0,1)^ and define r{x) := \x — x\, where x = (1/2,1/2). We choose r = 0.1, m = 0.35, 
M = 0.45 and cr = 1, and define the (optimal) state 


y{r) 


- 0 . 1 , r<l, 

< 614.4r^ - 768r^ + 352r3 - 72r^ + ^ ^ 

0 r > - 


and the adjoint state 


q{r) 
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- 0.1 
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0 0 
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0 0 



0 0 


0.02 

0.015 



0 0 


Figure 1: Optimal (top) and optimal discrete piecewise constant (bottom) states (left), 

controls /, (middle) and adjoint states g, (right) with 7 = 10^ and h — Hqui Example 1. 


The multiplier v is composed of a regular part concentrated in := S(x, 1/8), and a part 
concentrated on the boundary Taking 


yair) = 


-5.1, r<l, 

(^i/(r), r > 1, 

we obtain as action of z/ G applied to an element g G C(fi) 


[ qdiy^ — [ qdx—^[ qds. 

Jn Jqi 4 


The auxiliary state z and the corresponding load / are set to 

z{x) = sin( 7 rxi) sin( 7 rx 2 ), and 

f{x) = —Az{x) = 27r^ sin( 7 rxi) sin( 7 rx 2 ). 

The optimal control I is given by 

-3/4 


l{x) = (P[^ 4 ^M 4 ]( 3 (gor)(a;) 2 :(x))) 


and 


enix) := - A{y o r){x) - z{x)l{x). 

The control and state variables are depicted in Figure where we set 

hk = \/2-2^-^, k>l. 
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Figure 2: Errors in the states (left) and controls (right) plotted against 7 for different 
values of h in Example 1 with piecewise constant state approximation. 


In all our numerical tests we approximate integrals using a 3-point Gauh quadrature rule. 
In the case of piecewise linear elements, however, we distinguish between those parts of the 
triangles on which the involved projections are active and those parts where they are inactive. 


Let us begin with the mixed state approximation. Table contains L^(ri)-errors of the 
state and L^(f])-errors of the control variables for a run of the path-following Algorithm]^ 
using the coupling 7 = 0(/i“^), starting with 74 = 400 at /i = /i 4 . We observe erro rs of the 
size 0(/i), or equivalently 0 ( 7 “^/^), which is twice the rate predicted by Theorem |4.6[ Figure]^ 
shows the development of the errors in the state and the control variables, respectively, over a 
large range of regularization parameters. It can be seen from the graphs that the discretization 
error for both variables is approximately of order one, which explains the above convergence 
rate of size 0 (/i), and is the square of the expected error bound derived in the previous 
section. Furthermore, the regularization error appears to be of order 0 ( 7 “^) for the states 
and between 0 ( 7 “^'^) and 0 ( 7 “^'^) for the controls, while Theorem 3.3 predicts an order 
of 0 ( 7 “^-^^). Computation on a series of random unstructured meshes with grid sizes in the 
range of the considered uniform meshes results in the same convergence rates, and thus rules 
out superconvergence effects. The observed convergence rates might be explained by the high 
regularity of the constructed solution. 

With the path-following Algorithm in the appendix 3 to 4 Newton steps are needed to 
compute the numerical solution, where we use the tolerance 10“^. This result is achieved 
independent of the grid size of the underlying mesh and thus indicates mesh-independence of 
the algorithm. 

Using piecewise linear, continuous state approximations, however, we in this example are 
not able to sufficiently progress in the regularization parameter. Even for small values of 7 and 
with damped Newton steps the semismooth Newton iteration failed to converge. This may be 
due to large slopes contributed by the term eo and by the jump in y^. Similar observations are 
reported in [Gunther and Hinze 2011 ] . where an interior point solver is used to treat gradient 
constraints in elliptic optimal control. There a jump in the exact control leads to oscillations 
in the discrete approximation with piecewise linear, continuous finite elements, which results 
in convergence problems for the solver. 
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Table 1 : Errors and corresponding experimental orders of convergence EOC for the piecewise 
constant states and controls in Example 1. Parameters are coupled via 7 = 

74 = 400. 


hk 

8 

1 

EOCj; 

‘-‘1 

L2 

EOC/ 


3.45e-2 

- 

5.66e- 

-1 

- 

he 

1.74e-2 

0.99 

3.32e- 

-1 

0.77 

h-j 

8.80e-3 

0.98 

1.82e- 

-1 

0.87 

hg 

4.39e-3 

1.00 

9.47e- 

-2 

0.94 

hg 

2.19e-3 

1.00 

4.83e- 

-2 

0.97 




Figure 3: Errors in the states (left) and controls (right) plotted against 7 for different 
values of h in Example 2 with piecewise constant state approximation. 


Example 2: Here we consider problem (P) with parameters r = 0.01, m = 0.1 and M = 
0.2, so that the physical assumptions of a thin plate are satisfied. We again set Q = (0,1)^ 
and define the load 


/(xi,X2) : = 


-0.04, xi < ^ 

0 . 01 , xi > 


We consider the numerical solutions for h — diS reference solutions. 

Figure]^ shows the behaviors of the state and control errors in the piecewise constant case. 
The large magnitudes of the errors in the dual control I stem from the small bounds 

on the primary control u. In this example the regularization error in the control exhibits 
two consecutive convergence behaviors: Up to 7 = 10^ we observe the theoretically derived 
order 0 ( 7 “^-^^), which for larger values of 7 improves to 0 ( 7 “^-^). In this example the piece- 
wise linear, continuous state approximation works well with our path-following Algorithm 
In Table we report our numerical findings for the state and control variables with piecewise 
constant and piecewise linear and continuous state approximation, resp. The parameters are 


coupled according to Theorems 4.4 and 4.6, i.e., 7 = 0{h ^), starting with 7 = 400 on / 14 , 
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Table 2: Errors relative to the reference solution (^ref? ^ref) on grid hg and corresponding 
experimental orders of convergence for the states yj^ and controls in Example 2. Parameters 
are coupled via 7 = 74 = 400, and 7 = 0(/i“^), 74 = 16, in the mixed and the 

piecewise linear, continuous case, resp. 



p.w. constant ansatz 


p.w. linear ansatz 



hk 

2/ref yJi 

LOO 

kef - ll 

L2 


2/ref yJi 

ifl 

/ref - ll 

L2 

/l4 

1 . 22 e -2 

3.04e+2 

3.82e-2 

9.95e-^ 

) 

5.81e+2 

/l5 

6.30e-3 

2.24e+2 

1.42e-2 

3.42e-^ 

) 

3.61e+2 

Hq 

3.54e-3 

1.71e+2 

3.80e-3 

1 . 21 e-^ 

) 

2 . 01 e +2 

hr 

2.09e-3 

1.23e+2 

1.29e-3 

4.58e-c 

i 

9.81e+l 


0.96 


0.44 


1.43 

1.54 


0.69 



0.83 


0.39 


1.90 

1.50 


0.85 



0.76 


0.48 


1.56 

1.40 


1.03 



and 7 = 0(/i“^), starting with 7 = 16 on /i 4 . With this coupling, the errors in the controls 
are roughly of the predicted orders, whereas the state errors seem to converge at a faster rate. 
Large slopes now occur in the dual control /, due to the small control constraints on u and the 
asymptotically singular behaviour of the adjoint state near the state active set, which reduces 
to a point in the limit 7 ^ 00 . 

Figure displays the resulting approximations to the optimal state y and control u on 
the grid with h — hj. Similar to the experiment in |Arnantu et al. 2000 ( Figure 4-1] with 
changing sign of the load, we observe a bang-bang-like control with minimal thickness along 
the boundary. 

Let us comment on the active set shapes of the variationally discretized control variable /, 
compare the discussion in Subsection |4.1[ Since the approximation of the state with i?To- 
elements yields piecewise constant control approximations the boundary of the active set in this 
case follows the finite element mesh. For piecewise linear, continuous states the control active 
set is generally bounded by piecewise, non-degenerate hyperbolas. In Figure this is depicted 
for the numerical solutions of Example 2 with 7 = 500 on the mesh with h = h^. As expected, 
in the piecewise linear, continuous case the boundary of the upper active set is already well 
approximated on this coarse mesh, given the rather small regularization parameter. 


Appendix A: Semismooth Newton method 


In order to solve the subproblems (D^) and (D^ via a semismooth Newton method (cf., e.g., 
[Hintermiiller, Ito and Kunisch 2003] ), we rewrite the system of primal and adjoint equations 
of the associated optimality systems, (4.4)-(4.5) and (4.12)-(4.13), respectively, in the matrix 
form 




A ki 
A 


(x^) = 0, 


(A.l) 
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0 

- 0.002 

- 0.004 

- 0.006 

- 0.008 

- 0.01 

- 0.012 

1 



0 0 



1 



0 0 



Figure 4: Final states (top row) and corresponding controls (bottom row) for h = hj 
in the path-following runs of Example 2. Piecewise constant (left) and piecewise linear and 
continuous (right) state approximation for 7 = 2.5600e+4 and 7 = 6.5536e+4, resp. 



Figure 5: Comparison of the upper active set boundaries of the control I in Example 2 with h = 
/i 5 and 7 = 500 in the Raviart-Thomas (left) and the piecewise linear, continuous case (right). 
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where the variational inequality, in the equivalent form of the projection formula, has been 
substituted for the control in the primal equation. The matrix A denotes the according 
finite element system matrix of the respective discretization ansatz, whereas the operators ki 
and ^2 result from the right-hand sides of the primal and adjoint equations, respectively. 
Due to the underlying projections these operators will not be Frechet-differentiable. Similarly 
to [Gunther and Tber 2009| we define the following generalized derivative of as 




A Dki{x.) \ 

Dk]{^) ^ ) ’ 


with D denoting the generalized derivative of Clarke (cf. [Clarke 1983| ). 

In the Raviart-Thomas case we use from [Bahriawati and Carstensen 2005] the Matlab 
code EBmfem in order to discretize the Poisson equations and refer the reader to this refer¬ 
ence for further details. With m and n denoting the numbers of edges and elements of the 
triangulation, let then ^ ]^2mn combined primal and dual state 

vector, and x = {^h^Vh^^q.h^OihY' ^ another vector. We have 


where 




kliqlT) 


klivlT) 


Zh\T{‘^ (Ph\T ^Mt) , Tei{qJ^), 

< M-^\T\zh\T, Tel{ql), 

^m-^\T\zhiT, Teu{ql), 

f 7(yjT +'^) 1^1 ’ Tea{yl), 

[ 0, else. 


with the control and state active and inactive sets 


(A.2) 


i{ql) := {Te7h\m^< 3qlzh < } , 

•= {T e7h\ Sqlzh <m^} , 

l{ql) -.= {Te7h\M^< 3qlzh } , 


o-iVh) ■= {T e7h\yl + T <0} . 


The generalized derivatives of ki and k^ are given by the diagonal matrices 


where 


Dki{x) := dmg{Dki{qh,Tj)).^^^ 

Dkli^) := d\^g{Dkl{yh,T^))^^^^ 


Dki{qh,T) 

Dkl{yH,T) 


1^1 Tei{qh), 

0 , else. 


7|r|, Tea{yh), 

0 , else. 


(A.3) 
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Algorithm 1 Sketch of the path-following method for solving problem (D). 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


^0^70 ^ positive initial grid size and regularization parameter 

xo ^ given initial state vector of appropriate dimension, with subvectors 

n ^ 0 

loop 

Zn ^ (scalar component of) Gh^if) 
while ||F^’^(x^)|| > given tolerance do 

end while 

In ^ (^[m4,M4] (3^n^n)) 

if + 7n/2||(yn + '^)“||i 2 (Q) is Sufficiently small then 

Stop and return the last iterate 
else 


^n+l ^ ^n/2 

jn+i ^ 7n * 2^ (with K. chosen as suggested by Theorems 4.4 and 4.6) 
x^+i ^ interpolation of x^ on the refined mesh 
n ^ n -\- 1 

end if 
end loop 


In the case of piecewise linear elements, let N be the number of inner nodes in the tri¬ 
angulation and a canonical basis of Y^. With x^ = G and an arbi¬ 

trary X = (y/i, q^Y ^ R^^ we have that 


Hivl) ■= (~ [ ^ ^(-oo,0] {Vh + 't) 

\ Jft / 

Denoting by g\ G and ^2 G 9P(_oo,o] subgradients of the projection operators, we 

further obtain 


j=i,-W 


(A.4) 


f 9 

Dki{^)j^k-= / -zl{Pirni^M^]{3qhZh))~'^^'^ gi{SqhZh)((>k((>jdx, 

J fl 

D^2i^)j,k ■= - 'y92{yh + r)(t)k(t>jdx. 

JQ 


(A.5) 


The nonlinear equation (A.l) can now be solved with a semismooth Newton method. For 
this purpose, and an initial iterate xq, we generate a sequence of semismooth Newton steps 
via 


x„+i Xn - DF^{xn) ^F^{xn), n G IN, 

which, very similar to Proposition 4.1 of [Gunther and Tber 2009| . can be shown to be well- 
defined and locally convergent. 

In order to approximate the solution of (D) one can perform a path-following method, as 
outlined in AlgorithmTo this end, the above iteration can be continued until the Euclidean 
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norm of the vector F^(x^), which measures how much violates the optimality system, is 
below a given tolerance, and use the obtained approximate solution to (D^) as the starting 
iterate on a refined mesh with increased 7 . This nested iteration approach helps to stay 
within the convergence radius of the Moreau-Yosida-penalized Newton method. Overall, in 
our computations we observed a marked sensitivity of the Newton-type method with respect 
to the regularization parameter. Should the Newton solver diverge, however, we were able 
to return into the domain of convergence by choosing a more conservative increase in 7 . As 
mentioned before, the associated controls can then be extracted from the state variables x^ 
by means of the projection formula (4.7). 
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